Preparation

Packages

library(tidyverse) 
library(openxlsx) #read excel file
library(ggrepel) #label text
library(vegan) #Biodiversity Index
library(sf) #Map
library(scales)
library(treemapify)

Call custom function

source("Custom_function.R")

Read data

Use read.xlsx() function to load data from Excel file to R.

#library(openxlsx)
filename <- "Plot_NC-BD_20230925_v4.xlsx"
# Bidoup plot
BDplot <- read.xlsx(filename, sheet = "BD_census2")

This dataset is 2 hectares:

BD-1 (1ha) is concentrated (100x100m)
BD-2 (1ha) is taken from 25 scattered subplots (20x20m).

Community (Species composition, size-class structure)

Number of family, genus and species

Use n_distinct() to counts the number of distinct value in column: family, genus, scientificName. And group_by() to count at each plot.

BDplot |>  
  group_by(PlotID) |> 
  summarise("family" = n_distinct(family),
            "genus" = n_distinct(genus),
            "species" = n_distinct(scientificName),
            "stem" = n())

## # A tibble: 2 × 5
##   PlotID family genus species  stem
##   <chr>   <int> <int>   <int> <int>
## 1 BD-1       46    70     124 12207
## 2 BD-2       53    94     169 11575

Forest structure

We usually group trees into 3 size-class of DBH:

small tree [1-5)cm;
medium tree [5-20)cm;
big tree >=20 cm

BDplot <- BDplot |>
  mutate(
    #convert dbh from mm to cm
    dbh_cm = dbh / 10,
    gDBH = case_when(dbh_cm >= 20 ~ ">=20", 
                     dbh_cm >= 5 ~ "[5-20)", 
                     TRUE ~ "[1-5)")
  )

head(BDplot)

We count the number of stem at each diameter class and each plot

(stemByDBH <- BDplot |>
   group_by(PlotID, gDBH) |> 
   summarise(n = n()) |> 
    mutate(freq = round(n / sum(n), 3)) |>
    arrange(desc(n)))

## # A tibble: 6 × 4
## # Groups:   PlotID [2]
##   PlotID gDBH       n  freq
##   <chr>  <chr>  <int> <dbl>
## 1 BD-1   [1-5)   9736 0.798
## 2 BD-2   [1-5)   9478 0.819
## 3 BD-1   [5-20)  2057 0.169
## 4 BD-2   [5-20)  1749 0.151
## 5 BD-1   >=20     414 0.034
## 6 BD-2   >=20     348 0.03

Figure

ggplot(BDplot, aes(factor(gDBH), fill = factor(PlotID))) +
  geom_bar(position = "dodge") +
  theme_bw(14) +
  labs(x = "DBH (cm)", y = "Nb of stem", fill = "PlotID")

We can also count the number of species at each diameter class.

(
  SpByDBH <- BDplot  |>
    group_by(PlotID, gDBH) |>
    summarise(
      "family" = n_distinct(family),
      "genus" = n_distinct(genus),
      "species" = n_distinct(scientificName)
    )
)

## # A tibble: 6 × 5
## # Groups:   PlotID [2]
##   PlotID gDBH   family genus species
##   <chr>  <chr>   <int> <int>   <int>
## 1 BD-1   >=20       28    35      53
## 2 BD-1   [1-5)      45    69     120
## 3 BD-1   [5-20)     42    61      96
## 4 BD-2   >=20       33    45      64
## 5 BD-2   [1-5)      53    91     160
## 6 BD-2   [5-20)     48    74     121

Figure

Comparison of family numbers between two plots

SpByDBH |> 
  ggplot(aes(x = gDBH, y = family, fill = PlotID)) +
  geom_col(position = position_dodge()) +
  geom_text(aes(x = gDBH, y = family, label = family),
            position = position_dodge(width = 0.9)) +
  theme_bw(15) +
  labs(x = "DBH class", y = "Number", fill = "Plot")

Comparison of genus numbers between two plots

SpByDBH |> 
  ggplot(aes(x = gDBH, y = genus, fill = PlotID)) +
  geom_col(position = position_dodge()) +
  geom_text(aes(x = gDBH, y = genus, label = genus),
            position = position_dodge(width = 0.9)) +
  theme_bw(15) +
  labs(x = "DBH class", y = "Number", fill = "Plot")

Comparison of species numbers between two plots

SpByDBH |> 
  ggplot(aes(x = gDBH, y = species, fill = PlotID)) +
  geom_col(position = position_dodge()) +
  geom_text(aes(x = gDBH, y = species, label = species),
            position = position_dodge(width = 0.9)) +
  theme_bw(15) +
  labs(x = "DBH class", y = "Number", fill = "Plot")

Tree distribution map

BDplot.BD1 <- BDplot |> filter(PlotID == "BD-1")

Draw a blank plot with draw_blank_plot_20ha() function.

p <- draw_blank_plot_20ha()
p

Add add tree to map

p + geom_point(data = BDplot.BD1,
               aes(x = gx, y = gy, size = dbh_cm, color =  gDBH),
               alpha = 0.7) +
  xlim(80, 180) + ylim(400, 500) +
  labs(color = "DBH class", size = "DBH(cm)")

For big tree

BDplotA <- BDplot.BD1 |> filter(gDBH == ">=20")
p + geom_point(data = BDplotA, 
               aes(x = gx, y = gy, size = dbh_cm),
               alpha = 0.7, shape = 21,fill = "blue", color = "white") +
    xlim(80, 180) + ylim(400, 500) +
    labs(size = "DBH(cm)")

For medium tree

BDplotB <- BDplot.BD1 |> filter(gDBH == "[5-20)")
p + geom_point(data = BDplotB, 
               aes(x = gx, y = gy, size = dbh_cm),
               alpha = 0.7, shape = 21,fill = "blue", color = "white") +
    xlim(80, 180) + ylim(400, 500) +
    labs(size = "DBH(cm)")

For small tree

BDplotC <- BDplot.BD1 |> filter(gDBH == "[1-5)")
p + geom_point(data = BDplotC, 
               aes(x = gx, y = gy, size = dbh_cm),
               alpha = 0.7, shape = 21,fill = "blue", color = "white") +
    xlim(80, 180) + ylim(400, 500) +
    labs(size = "DBH(cm)")

For a particular species, px: Pinus krempfii

## So do phan bo toan bo cay trong o mau
BDplotSP <- BDplot.BD1 |> filter(scientificName == "Pinus krempfii")

p + geom_point(data = BDplotSP,
               aes(x = gx, y = gy, size = dbh_cm, color =  gDBH),
               alpha = 0.7) +
  xlim(80, 180) + ylim(400, 500) +
  labs(color = "DBH class", size = "DBH(cm)")

We can also try with a detail DBH class

BDplot.BD1 <- BDplot.BD1 |> 
mutate(
  dbh_cm = dbh / 10 , # dbh in mm
  gDBH2 = case_when(
                dbh_cm >= 100 ~ ">=100",
                dbh_cm >= 90 ~ "[90-100)",
                dbh_cm >= 80 ~ "[80-90)",
                dbh_cm >= 70 ~ "[70-80)",
                dbh_cm >= 60 ~ "[60-70)",
                dbh_cm >= 50 ~ "[50-60)",
                dbh_cm >= 40 ~ "[40-50)",
                dbh_cm >= 30 ~ "[30-40)",
                dbh_cm >= 20 ~ "[20-30)",
                dbh_cm >= 10 ~ "[10-20)",
                TRUE ~ "[1-10)")
    )
# Count number of stem in each DBH class
stemByDBH2 <- BDplot.BD1 |>
  count(gDBH2) |>
  mutate(freq = round(n / sum(n), 3)) |>
  arrange(desc(n))

stemByDBH2

##       gDBH2     n  freq
## 1    [1-10) 11169 0.915
## 2   [10-20)   624 0.051
## 3   [20-30)   246 0.020
## 4   [30-40)   123 0.010
## 5   [40-50)    35 0.003
## 6   [70-80)     3 0.000
## 7     >=100     2 0.000
## 8   [50-60)     2 0.000
## 9   [60-70)     1 0.000
## 10  [80-90)     1 0.000
## 11 [90-100)     1 0.000

# Figure
stemByDBH2 |>
  ggplot(aes(x = gDBH2, y = n)) +
  geom_col() +
  geom_text(aes(label = n), vjust = -0.5) +
  theme_bw(14) +
  labs(x = "DBHClass (cm)", y = "Number of stem")

Dominant species and Biodiversity Index

Dominant species

The Importance Value Index (IVI) in Ecology is the quantitative measure of how dominant a species is in a given ecosystem. IVI is calculated from 3 values:

relative density of each species (RD)
relative basal area of each species (RBA)
relative frequency of each species (this can be omitted in 1 plot)

\[IVI={\left(relative.density + relative.basal.area \right)}\]

Basal area (BA) of each species can be calculated from diameter measured: $BA = pi * R^2$

Relative basal area (RBA): $RBA = BA / sum(BA)$

Relative density (RD) = number of individuals of each species / total

  BDplot <- BDplot |> 
              mutate(BA = ((dbh/2)/1000)^2*pi) # unit: m2)

BDplot.BD1 <- BDplot |> filter(PlotID == "BD-1")
BDplot.BD2 <- BDplot |> filter(PlotID == "BD-2")

IVI.1 <- calcIVI(BDplot.BD1)
head(IVI.1)

## # A tibble: 6 × 7
##   scientificName              n    sG     UT    MD   IVI  rank
##   <chr>                   <int> <dbl>  <dbl> <dbl> <dbl> <dbl>
## 1 Lithocarpus pierrei       966 6.74  12.7   7.91  20.6      1
## 2 Polyosma nhatrangensis   1175 2.33   4.40  9.63  14.0      2
## 3 Syzygium rubicundum       726 3.21   6.06  5.95  12.0      3
## 4 Pinus krempfii             89 5.28   9.96  0.729 10.7      4
## 5 Castanopsis echinocarpa   539 3.05   5.75  4.42  10.2      5
## 6 Lasianthus saprosmoides  1076 0.236  0.444 8.81   9.26     6

IVI.2 <- calcIVI(BDplot.BD2)
head(IVI.2)

## # A tibble: 6 × 7
##   scientificName              n    sG    UT    MD   IVI  rank
##   <chr>                   <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Lithocarpus coinhensis    357 3.94  7.84   3.08 10.9      1
## 2 Castanopsis echinocarpa   483 3.27  6.51   4.17 10.7      2
## 3 Polyosma nhatrangensis    705 1.31  2.60   6.09  8.69     3
## 4 notdet                    883 0.127 0.253  7.63  7.88     4
## 5 Rhaphiolepis poilanei     436 1.07  2.13   3.77  5.90     5
## 6 Lithocarpus pierrei       262 1.40  2.79   2.26  5.05     6

Figure

figureIVI(IVI.1)

figureIVI(IVI.2)

Diversity Index

biodivIndex(IVI.1)

##   Richness   Simpson  Shannon  Evenness
## 1      124 0.9598452 3.765398 0.7811574

biodivIndex(IVI.2)

##   Richness   Simpson  Shannon  Evenness
## 1      169 0.9756952 4.200136 0.8187561

Biomass and carbon stock

Biomass

Above-ground biomass (AGB)
Below-ground biomass (BGB)

AGB calculation

Chave et al (2005). Tree allometry and improved estimation of carbon stocks and balance in tropical forests. Oecologia, 145(1): 87–99.

AGB=Total aboveground biomass (kg)
D=DBH (cm)
H=Height (m)
$ρ$=Wood density (g/cm3)
BA=Basal area (cm2)

Habitat	Equation in Chave et al (2005)
Tropical forest DRY with Height	$AGB = exp(-2.187 + 0.916 * ln(ρ * D^2 * H)$
Tropical forest DRY only DBH	$AGB = ρ * exp(-0.667 + 1.784 * ln(D) + 0.207 * (ln(D))^2 – 0.0281 * (ln(D))^3)$
Tropical forest MOIST with Height	$AGB = exp(-2.977 + ln(ρ * D^2 * H))$
Tropical forest MOIST only DBH	$AGB = ρ * exp(-1.499 + 2.148 * ln(D) + 0.207 * (ln(D))^2 – 0.0281 * (ln(D))^3)$

library(BIOMASS)  

woodensity <- getWoodDensity(
  genus = BDplot.BD1$genus,
  species = BDplot.BD1$scientificName,
  family = BDplot.BD1$family
)

head(woodensity)

BDplot.BD1$meanWD <- woodensity$meanWD

Bidoup Plot: Tropical forest MOIST only DBH

$AGB = ρ * exp(-1.499 + 2.148 * ln(D) + 0.207 * (ln(D))^2 – 0.0281 * (ln(D))^3)$

BDplot.BD1 <- BDplot.BD1 %>% 
  mutate(dbh_cm = dbh/10,
         AGB = meanWD * exp(-1.499 + 2.148*log(dbh_cm) + 0.207*log(dbh_cm)^2 - 0.0281*log(dbh_cm)^3))

(AGB <- sum(BDplot.BD1$AGB)/1000) #Ton/ha

## [1] 469.1828

BGB calculation

The ratio developed by Mokany et al. (2006). Critical analysis of root : shoot ratios in terrestrial biomes. Global Change Biology, 12: 84-96 offers specific ratios based on forest type and climate zone and are applicable when the aboveground biomass estimate (shoot) is reported at the stand level

(BGB = 0.235 * AGB)

## [1] 110.2579

Total biomass of 1 hecta

(TB <- AGB + BGB)

## [1] 579.4407

Carbon stock

The quantity of carbon can then be estimated by converting biomass to carbon using the IPCC default carbon fraction of 0.47. After that, the carbon stock was multiplied by 44/12 (the carbon atom ratio in the molecular weight of CO2) to get the CO2 equivalent value.

#Carbon stock (ton/ha)
(TC <- TB * 0.47)

## [1] 272.3371

# CO2 (ton)
(CO2 <- TC * 44/22)

## [1] 544.6743

Carbon credit

1 carbon credit = 1 ton of CO2 equivalent (tCO2e) emissions.

In 2023, Vietnam successfully sold 10.3 million forest carbon credits for the first time through the World Bank at a price of $5/ton (Global average selling price ~$11.2)

1 hecta forest at Bidoup - Nui Ba national park ~ 2700 $

(CO2 * 5)

## [1] 2723.371

Tree distribution map in large scale

LITHCA : Lithocarpus eucalyptifolius
CASTEC: Castanopsis echinocarpa
SYMPBA: Symplocos banaensis

Day 02 - Bidoup Forest Dynamics Plots

Nguyen-Le Xuan Bach

2025-11-30

Preparation

Packages

Call custom function

Read data

Community (Species composition, size-class structure)

Number of family, genus and species

Forest structure

Tree distribution map

Dominant species and Biodiversity Index

Dominant species

Diversity Index

Biomass and carbon stock

Biomass

AGB calculation

BGB calculation

Carbon stock

Carbon credit

Tree distribution map in large scale